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Abstract 

We propose a way of implementing the Broad Histogram Monte 
Carlo method to systems with continuous degrees of freedom, and we 
apply these ideas to investigate the three-dimensional XY-model with 
periodic boundary conditions. We have found an excellent agreement 
between our method and traditional Metropolis results for the energy, 
the magnetization, the specific heat and the magnetic susceptibility on 
a very large temperature range. For the calculation of these quantities 
in the temperature range 0.7 < T < 4.7 our method took less CPU 
time than the Metropolis simulations for 16 temperature points in that 
temperature range. Furthermore, it calculates the whole temperature 
range 1.2 < T < 4.7 using only 2.2 times more computer effort than 
the Histogram Monte Carlo method for the range 2.1 < T < 2.2. Our 
way of treatment is general, it can also be applied to other systems 
with continuous degrees of freedom. 



1 Introduction 



For thermodynamic systems in equilibrium one typically wants to calculate 
average values of quantities Q like the energy, the magnetization, etc.. These 
mean values depend on one hand of the properties of the system himself and 
on the other hand on its interactions with the environment. For example, in 
the case of the canonical ensemble, the system is in equilibrium with a heat 
bath at temperature T, and the mean value < Q >t is given by 

j:E9{E)exp{-E/ksT) ' 

where the degeneracy g{E) is the number of distinct states with energy E, 
< Q >E denotes the micro-canonical average of Q at energy E and ks is the 
Boltzmann constant (in the rest of this paper fc^ = 1 is taken). Here, g{E) 
and < Q >e are intrinsic characteristics of the system, and exp {—E/kBT) 
is the Boltzmann factor. 

Most Monte Carlo methods were designed to calculate particular averages 
like Eq. (0). For example. Metropolis algorithms and similar canonical 
simulations accumulate sampling distributions that approximate Eq. (|l]) for 
a large number of samples and, therefore, the mean values < Q >t are 
estimated simply by averaging Q on the sample set. However, if one wants 
to study the dependence of < Q >t on T using this kind of algorithms it is 
necessary to run one entire simulation at each temperature value, with large 
computer efforts. 

A well-proved strategy to avoid these multiple calculations is the his- 
togram method (HMC), introduced by Salsburg p[ and popularized by Fer- 
renberg and Swendsen [Q. This method reweights data of one canonical 
simulation at one temperature To to find mean values at other temperatures. 
However, a canonical distribution has a very narrow peak on the energy axis 
around its mean energy and, therefore, there are not enough samples in the 
tails of the distribution to have good statistics there. This restricts the re- 
liability of the method to a relatively narrow temperature range around To 

The Broad Histogram Monte Carlo Method (BHMC) ^-[1^ was estab- 



lished by de Oliveira et. al. in order to overcome these restrictions. It is 
designed to determine the system density of states g{E) and < Q >e di- 
rectly. These quantities are calculated from micro-canonical averages and. 
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therefore, samples with different energies can be considered independently. 
The samples can be obtained in many different ways, for instance, perform- 
ing micro-canonical simulations |^ or performing a non-biased random walk 
along the energy axis, which gives rise to much broader histograms than the 
usual histogram method P, |^. Given g{E) and < Q >e, the canonical dis- 
tribution and the desired mean values for any temperature can be obtained. 
Thus, it is possible to calculate mean values over the whole temperature 
range of interest in only one run. The BHMC method was first developed for 
systems with discrete degrees of freedom, and has shown to give results more 
accurately and using less computer time than the conventional histogram 
method for the 2d an 3d Ising model the Edwards Anderson spin 

glass [0 and the 3d Potts glass 

In the present paper we propose a way of extending the method to systems 
with continuous degrees of freedom, using the XY-model as testing ground. 
Our way of treatment is general, it can also be applied to other systems 
of this kind. The remaining part of the paper is organized as follows. In 
Sec. we briefly summarize the basis of the BHMC method. In Sec. |^ we 
describe the XY-model and we use it as example to show how to apply the 
method to continuous systems. In Sec. ^ we discuss several ways of taking 
the samples. In Sec. ^ we give technical details about the implementation we 
did of the method. In Sec. ^we present the results obtained with the BHMC 
method from a 3D XY-model on a 10 x 10 x 10 cubic lattice with periodic 
boundary conditions, and we compare these results with those of Metropolis 
simulations and of the histogram method. On one side, we have obtained an 
excellent agreement between the BHMC method and Metropolis simulations 
in all the calculated quantities over a very broad temperature range, i.e. for 
temperatures between 0.7 and 4.7, using less computer time for the BHMC 
method than the time required for 16 points by Metropolis simulations. On 
the other side, it reproduces with the same accuracy the results of the HMC 
method on the temperature range 2.1 < T < 2.2 but, instead of diverging for 
temperatures out of this range as the HMC method does, it remains precise 
over the whole range 1.2 < T < 4.7 using only 2.2 times more computer 
effort than the HMC method. It justifies the name "Broad Histogram" for 
the method. Finally, in Sec. ^ we summarize the main steps of the proposed 
strategy to apply the BHMC method to continuous systems and we discuss 
other strategies that can be used to implement the method. 
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2 The BHMC Method 



The strategy of the BHMC method is to calculate the degeneracy of energy 
states g{E) directly. The method for systems with discrete degrees of freedom 
can be summarized as follows. 

First, let us imagine a protocol of allowed movements in the space of states 
of the system such that changing from an Xgid to an X„e«, configuration is 
allowed if and only if the reverse change is also allowed, i.e. 

-^old -^new -^new ^ -^old (2) 

is allowed is allowed 

In other words, the protocol is micro-reversible. These movements are only 
virtual, in the sense that they are not performed. They are introduced only 
to estimate the density of states g(E). 

Next, let us choose a fixed amount of energy change AEfi^ and com- 
pute Nup{X) (iVrf„(X)) for the configuration X as the number of allowed 
changes that increases (decreases) the energy of the configuration by AEfi^. 
Let < Nup{E) > (< Ndn{E) >) be the micro-canonical average of Nup{X) 
(NaniX)) at energy E. 

Due to Eq. (^, the total number of ways to go up, summed over all states 
with energy E, is equal to the total number of ways to go down, summed 
over all states with energy E + AEfi^, i.e. 

g{E) < N^p{E) >= g{E + AEf,^) < Ndn{E + AE/,,) > . (3) 

Eq. d^) can be rewritten as: 

Hence, knowing < Nup{E) > and < N^niE) > allows to calculate the right 
hand side of Eq. (^) and obtains In g{E) for all values of E in steps of AEfi^ 
by adding up differences from a given initial point {Eo,ln g{Eo)). 

Similarly, it can be observed that Eq. (H) divided by AEfi^ estimates 
f3{E), i.e. 

dlng{E) lng{E + AEfi,)- In g{E) 



(3{E) 



dE AE 



fix 



1 |„ <N,AE)> 



AEf,^ < Ndn{E + AEf,.,)> 
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Then, P{E) can be obtained for a set of energy values {Ei, E2, ■ ■ ■ , En} by 
estimating < Nup{Ei) > and < Ndn{E + AEi f^^^) > for each value Ei. Eq. (^) 
can be used to estimate P{E) in more accurate ways, as it will be discussed in 
Sec. |[ \ng{E) can be obtained, thereafter, by integrating P{E) numerically 
from a given initial point {Eo,ln g{Eo)). 

The value of In g{Eo) for the initial point does not affect the calculation 
of mean values < Q >t, as can be shown in Eq. (|l|). If In g{Eo) = is taken, 
the method gives g{E)/ g{Eo) for each energy. 

Finally, to carry out calculations with this method four histograms are 
required, i.e. Nup{E), Ndn{E), Q{E) and the number of visits V{E). Estima- 
tions for < Nup{E) > and < Ndn{E) > and, in consequence, an estimation 
of g{E) are obtained by dividing Nup{E) and Ndn{E) by V{E). < Q >e is 
calculated by dividing Q{E) by V{E). Combining these two results Eq. ([^) 
can be used to calculate < Q >t at any desired temperature. 

It can be shown that Eq. (^) is exact for any statistical model |T0|. Monte 



Carlo sampling processes are used only to estimate the micro-canonical aver- 
ages < Nup{E) > and < Ndn{E) >. Many different sampling strategies can 
be employed, as will be discussed in Sees. | and ^. Nevertheless, all states 
with the same energy have to be sampled with the same probability and the 
correlations between samples have to be small enough to obtain good esti- 
mations of such micro-canonical averages. Otherwise, problems can appear 
0,11011. 



3 The BHMC method and the XY-Model 
3.1 The XY-model 

In order to describe our technique for extending the BHMC method to sys- 
tems with continuous degrees of freedom we have chosen the XY-model as 



testing ground. The XY-model |T2[-|1^, consists of a set of spins a of length 



unity arranged on a lattice. Each spin is allowed to rotate in a plane, charac- 
terized by an angle 6. These variables define the state of the system and can 
change continuously in the range 6 G [— vr, tt]. The Hamiltonian with zero 
external field is given by 

n = -JY.a,-a, = -jY.cosie,-ej) (6) 

<ij> <ij> 
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where the summations < ij > are taken over all pairs of nearest-neighbor 
sites, (Tj ■ cTj is the scalar product between a, and aj, and J denotes the 
maximal energy per bond (in the rest of this paper we will take J = 1). The 
order parameter chosen for this model [|l^] is the average magnetization M, 
defined by: 

M^^A, , (7) 



with 



\ 



Jocose J + 51 sin 



(8) 



where the summations j are taken over all nearest neighbors of the site i. 
The XY-model so defined is one of the simplest thermodynamic systems with 
continuous degrees of freedom. 



3.2 The BHMC method for continuous systems 

In order to extend the BHMC method for continuous systems, it is necessary 
to redefine the quantities involved in Eq. and to find a condition for the 
protocol of movements equivalent to Eq. (||) such that Eq. (j^) holds. 

First, let us imagine a protocol of random movements in the space of 
states of the system, such that for each allowed movement the probability 
to perform it equals the probability to revert it. These movements are again 
only virtual, in the sense that they are not executed. They are introduced 
only to estimate the density of states g(E). For the XY-model, for example, 
we use the following protocol: we choose one spin at random with angle doid, 
and then we choose for it a new angle 6new according to a uniform probability 
distribution over the range 9new £ [— ^r, +7r]. It has to be noted that 9new is a 
random variable, with a probability density function (p.d.f) fe„^^,{^) defined 
as the probability of obtaining a value 6new between 6 and 6 + d6 given by 

J0ne^ If^J - j Q . otherwise ' 

Thus, the protocol of random movements can be defined by giving the p.d.f. 's 
for the new values of the system variables. This seems a natural implemen- 
tation for continuous systems. 
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The next step consist in estimating Nup and Ndn for a given configuration 
with energy Eoid- The problem arises because there is a continuum non- 
numerable set of configurations that can be reached. However, it has to be 
noted that each one of these new configurations has a well-defined energy 
value Enew and that exactly one of them is chosen at random. Therefore, 
we can define Enew and the energy change AE = Emw — Eoid as random 
variables. The p.d.f. fAE^AE) is defined by the probability of obtaining an 
energy change between AE and AE + dAE. We propose to redefine N^p and 
Ndn as 

Nup = fAE{AEf,,) ; Ndn = fAEi-AEfi,) (10) 

From this point, the method proceeds the same as in the discrete case. Eqs. 
( p!oD extend the BHMC method to systems with continuous degrees of free- 
dom. Therefore, the problem reduces to finding the function /ae for a given 
configuration Xoid- 

3.3 Finding f^E 

The procedure of finding fj\E is straightforward. First, the protocol of ran- 
dom movements is defined by giving the p.d.f. 's for the new values of the 
system variables, as in Sec. Then, the energy Emw, is expressed in 

terms of such new values. Next, its p.d.f. fE^em is determined by using the 
usual rules of finding the p.d.f. for a function of random variables Fi- 
nally fAEiAE) is found, according to these rules, replacing E in /£;„<.^(-E') 
by AE + Em. 

To continue with our example, we analyze first what happens with the 
turn of one specific spin. Suppose that we choose to turn spin i. Its bonds 
are the only ones that change the energy value. We define Si as sum of the 
energies of these bonds, i.e. 

e^ = -Y^cos{e^-e,) . (11) 

where the summations are taken over all nearest neighbors j of the site i. 
When 6i changes from its previous value 6ioid to the new value Oinewi 
changes from ei^id to Sinew It is clear that AE = Sinew ~ ^ioid and that Sinew 
is a function of the random variable Oinew- Therefore, Si^id and Sinew can be 
used instead of Enew and Eoid for the general treatment described above. 
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In order to express Einew iii terms of 6'j„g^, Eq. ( ]Tl| ) can be rewritten as 
where 

bQi = arctan ^ + vr (13) 

22 j cos Oj 

and Ai is defined by Eq. (|]). 

If 6'j„e«) is a random variable with p.d.f. given by Eq. (H), the usual rules 
of finding the p.d.f. for a function of random variables can be used to show 



that Einew is a random variable with p.d.f. given by |]T5 

(e) = I ■ ^ . (14) 

[ : otherwise 

Finally, we replace e in Eq. (|1^ by /S.E + Sioid, to obtain the p.d.f. 
f^^lAE) of Ai^^ corresponding to a turn of the spin i, i. e. 

fU^E) = \ VA-feo.+Ai.p -'^1 . (15) 

I : otherwise 

Since all spins are equally probable, the probability of obtaining an energy 
change AE between AE and AE + dAE is equal to the sum of f}^^{AE) 
over all spins divided by A^, the number of spins, i.e. 

UE{AE) = j^Y.rAEi^E) . (16) 

i 

Nup and N^n for a given configuration are obtained finally using Eqs. (|l^) as 
Nup = l^EKp ; N,^ = j^Y.Nk (17) 

i i 

where we have defined 

K,^ flEi^Ef,,) ■ Nl^^flj^i-AEf,,) (18) 

When Eqs. ([17|) are replaced in Eq.(^, the factor in Eq. ([T6| ) as well 
as the factor l/vr in Eq. (|T5[) will eventually cancel out. 
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Figure 1: Probability Density function f\^{AE) of the energy change AE obtained by 
turning spin i, according to Eq. (|l^). The function diverges at points AE = ztAi — Eidd, 
takes its minimum value between these points at AE — —eioid and equals zero outside 
this range. The contributions A^^^ and A^^^^ of turning the spin i are shown. 



The p.d.f. f^^{AE) is shown in Fig. |l|. The function diverges at points 
AE = ±Ai — Eigid and takes its minimum value between these two points at 
AE = Sold with fAEi~^ioid) = (7rv4j)^^. Outside this range the function 
equals zero. The values N^p and N^,^ are also shown. 

Summarizing, Eqs. ( [T7| ) give us the way to calculate Nup and Ndn for the 
XY-model and, therefore, a way to apply the BHMC method to it. The same 
strategy can be used to apply the method to other continuous systems. 

4 How to take the samples 

A second idea behind the BHMC method is to take the samples in such a 
way that they distribute almost homogeneous over the whole energy range 
of interest. There are many different strategies that can be used for this 
purpose. We have used two of them, namely, a micro-canonical sampling 
process that maintains the system inside a narrow energy window (BHMC- 
/iC) [§], and a random walk on the energy axis using Metropolis ^ steps 
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(BHMC-M) §. 

The BHMC-/xC sampling strategy consist in taking a small window cen- 
tered at each energy of interest, starting from a configuration inside the 
window and performing random movements never surpassing the limits of 
the window. It can be implemented, for instance, by taking at random a 
new configuration and accepting it only if it falls inside the window, or by 
using the Creutz energy bag method fl^. According to Eq. (|), two micro- 
canonical simulations are required to calculate P{E) at one energy value. 
However, the data of one simulation can be used for two energy values if 
their difference equals AEfi^. For the XY- model, we take one spin and gen- 
erate at random a new angle 6 G [— 7r,7r]. If the total energy of the system 
falls inside the window, the change is accepted, otherwise it is rejected. It is 
clear that this sampling process maintains a detailed balance condition with 
equal probabilities for the old an new configurations and, therefore, it sam- 
ples all states inside the window with the same probability. A whole lattice 
sweep is performed by repeating this procedure for all spins of the system, 
and a new configuration is taken after every fixed number of lattice sweeps. 

The BHMC-M sampling strategy consist first in performing a Markovian 
process with symmetric probability distributions for increasing and decreas- 
ing the energy, i.e. a non biased random walk on the energy axis. Sec- 
ond, since the canonical distribution is almost symmetric on the energy axis 
around its mean energy value < E >t, one way to do it (but not the only 
one) consists of taking the last configuration of the sample and performing 
on it a fixed number of Metropolis [|l| steps at such a temperature that its 
canonical distribution appears centered on the energy value E. The resulting 
configuration is taken as the next sample. In order to find this temperature 
T{E) Eq. d^) can be used by remembering that l3{E) = 1/T[E). Further- 
more, Eq. (^) can be used to estimate (3{E) and T{E) in more accurate ways. 
In the present paper we use the centered difference estimator 

\ng{E + AEj,,) - \ngiE - AEj,,) 
P\E) ~ 



2AE 



fix 



, < Nup{E) >< Nup{E - AEfi,) > 

^ AT ^ ^ AT riP . /KTP y^^) 



2AEfi, < Ndn{E) >< Nan{E + AEf^,) > ' 

These values are estimated from the data accumulated in the histograms 
Nup{E), Ndn{E) and V{E). It can happens that the system reaches an 
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energy value E such that some of the neighboring energies have yet not 
been sampled and, therefore, T{E) cannot be estimated. In such the 
Metropolis steps are performed at the same temperature used to produce the 
last configuration. 

5 Implementation 

In order to test our strategy numerically we performed computations for the 
three-dimensional XY-model on a 10 x 10 x 10 cubic lattice with periodic 
boundary conditions using the BHMC-M, the BHMC-/xC, the HMC and the 
Metropolis method. 

For the Metropolis method the spin system was initialized at each tem- 
perature value from a random configuration performing iVj^j = 500 entire 
Metropolis lattice sweeps before sampling, and then A^^ = 1500 samples were 
taken separated by Ni = 10 lattice sweeps among them in order to decrease 
their correlations. Here and in the following the correlation between succes- 
sive samples was determined by looking at the correlation of their magneti- 
zation values. For the Metropolis method we obtained a correlation between 
successive samples of 56.5% at To = 2.159, i.e. at the pseudo critical temper- 
ature we have observed for this finite system. Some additional Metropolis 
simulations were performed as reference data for the comparison between the 
BIIMC-/iC and the HMC methods. These simulations were performed in the 
same way as before, but taking A^^ = 10000 samples at each temperature 
value. 

For the BHMC-/iC method we divided the whole energy range of posi- 
tive temperatures, i.e. energies per bond between —1.0 and 0.0, in adjacent 
windows of equal size and chose AEfi^ to be equal to this bin size. So 
AEfix is equal to the difference of the energy values at which consecutive 
micro-canonical simulations are performed and, therefore, the data of each 
simulation can be used twice (see Sec. To acquire the samples the spin 
system was initialized for the lowest-energy window from a random config- 
uration by turning some spins at random until the window is reached and 
performing Nini micro-canonical lattice sweeps before sampling. Then, we 
took Ng samples separated by Ni lattice sweep among them. After that, we 
move to the next window by turning some spins at random starting from the 
previous sample. Then, Nini micro-canonical lattice sweeps were performed 
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before sampling, and so on. 

To compare with the Metropohs simulations, we divided the whole energy 
range of positive temperatures in 1225 adjacent windows, i.e. AEji^ ~ 2.45. 
To handle the critical slowing down in the relaxation time for the micro- 
canonical simulations, larger values of N^ni have been taken in the vicinity 
of the critical point, i.e. Nini = 250 for energies per bond between —0.400 
and —0.320 and Nini = 70 for energies outside this range. To obtain the 
same correlation between samples as from Metropolis simulations we used 
A^j = 3 in the energy range (—0.429, —0.282) and Ni = 2 outside this range. 
The maximal correlation in the vicinity of the critical point so obtained was 
56.4%. We took A^^ = 50 for the whole energy range. 

To compare with the HMC method, the whole energy range of positive 
temperatures was divided in 500 adjacent windows, i.e. AEfi^ = 6.0. To 
handle the critical slowing down, we took Nini = 200 for energies per bond 
between —0.400 and —0.320, and Nini = 80 for energies outside this range. 
Ni = 2 was enough to obtain lower correlations in the vicinity of the critical 
point than from Metropolis simulations. The maximal correlation so obtained 
was 46.3%. In order to obtain the same accuracy than the HMC method on 
the temperature range 2.1 < T < 2.2 more samples have been taken in the 
vicinity of the critical point, i.e. Ns = 380 for energies per bond between 
—0.504 and —0.248 and Ns = 80 for energies outside this range. 

For the BHMC-M method the histograms were constructed by divid- 
ing the whole energy range of positive temperatures in 1225 boxes of equal 
size. Ten independent systems that accumulated data in common histograms 
were simulated simultaneously, in order to estimate T{E) more accurately. 
These systems were initialized from random configurations thermalized by 
Nini = 500 Metropohs lattice sweeps at the critical temperature Tc = 2.20196 



1^. We took a total number of iV^ = 10 x 10000 samples performing Ni = 1 
lattice sweeps between them. In order to avoid wasting samples at very low 
temperatures, we restricted our systems to energies between —0.90 and 0.0 
per bond; spin systems going out of this range were reinitialized by copying 
the values of an additional spin system, which acts as a replacement. This 
additional system was initialized in the same way as the other ones, and then 
twenty entire Metropolis lattice sweeps at were applied on it after each 
time it was taken as replacement. 

The HMC method was implemented by dividing the whole energy range of 
positive temperatures in 1225 and 500 boxes of equal size, and accumulating 
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the data of a Metropolis simulation at To = 2.159 performed with Nini = 
500,Ni = 10 and iV, = 10000. 

For all methods, we stored cos(^j) and sin(^j) in tables instead of the 
angles 6i in order to accelerate the calculations. For the BHMC-/iC and 
BHMC-M methods, we stored in addition and for each spin for the 
same reason. The functions cos(a;), sin(a;), a/x, l/^x, exp(x), were stored in 
tables for all methods. 

The mean values and error bars for each method were made by repeating 
eight times the whole simulation with different initial seeds of the random 
number generator. We used the "minimal" Park and Miller random number 



generator combined with a Marsaglia shift sequence from Ref. [17 



6 Results 

The results are shown in Fig. p|-|To|. Figs. display data obtained by using 
1225 divisions of the energy range of positive temperatures, i.e. AEfi^ ~ 2.45 
for the BHMC-/iC, the BHMC-M and the HMC methods. 

Fig. ^ shows how samples distribute on the energy axis using the BHMC- 
fiC, the BHMC-M and the HMC methods. For the BHMC methods a large 
energy range is sampled, in contrast to the narrow window sampled by the 
classical canonical histogram method. In addition, with the BHMC-/iC it is 
possible to take exactly the same number of samples at each energy value. 

Fig. I shows P{E) obtained from the BHMC-/iC and the BHMC-M meth- 
ods. It has to be noted that both curves are almost indistinguishable, in spite 
of the completely different sampling strategies of both methods. 

Fig. I shows \ng{E) obtained from the BHMC-^uC and the BHMC-M 
methods. Since the curve from the BHMC-M method was obtained adding 
the differences of Eq. (^) starting from the point {Eo = —0.90, \iag{Eo) = 0), 
its values g{E) are the ratio of the probability to obtain at random an energy 
per site between E and E + dE to the probability to obtain an energy per site 
between this value of Eq and Eo + dE, when all states are equally probable. 
The same applies to the curve obtained from the BHMC-/iC method, with 
Eo = —0.94. It can be noted, that both curves are only shifted in a constant 
value, namely ln(c/(-0.90)/^(-0.94)). 

Figs. I^-P compare results from the BHMC methods with Metropolis simu- 
lations. Figs. ^ and ^ show the average energy and the average magnetization 
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Figure 2: Number of visits as a function of energy obtained from the BHMC-/iC method 
with 50 samples per window, the BHMC-M method with 10 x 10000 samples, and the 
canonical histogram method with 10000 samples for the 10 x 10 x 10 XY-model on a cubic 
lattice. All three calculations were performed by using 1225 divisions on the energy axis. 
The points over and below each curve correspond to the upper and lower limits of the 
error bars. On this energy scale the mean energy at the critical temperature observed 
corresponds to E = —0.377. 



obtained with the BHMC-yuC and the BHMC-M methods compared with the 
same quantities obtained using Metropohs simulations. The agreement be- 
tween the three Monte Carlo methods is excellent. The data obtained by the 
BHMC methods fit so well with those obtained by Metropolis simulations, 
that their curves and their error bars are almost indistinguishable. 

In Fig. ^ and | we show the heat capacity and the magnetic susceptibility 
obtained by using the BHMC-/iC, the BHMC-M and the Metropolis simu- 
lations. The curves calculated by the BHMC-M method fit very well within 
their error bars with them of the Metropolis simulations for the temperature 
range between T = 4.7 and T = 0.7, and deviate for lower temperatures. 
The curves calculated by the BHMC-/iC method fit as well as those from 
the BHMC-M method in the same temperature range, begins to oscillate at 
T = 0.7 and deviates for a lower temperature, namely T = 0.45. This is due 
to a lack of good statistics at low energies. In the case of BHMC-/iC method 
most of the possible changes at low energy values increase the system energy 
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Figure 3: /3 = l/T as a function of the energy for the SC 10 x 10 x 10 XY- model obtained 
from the BHMC-/.tC method (dotted hne) and the BHMC-M method (sohd Hne). Both 
calculations were performed with AEfix — 2.45. The error bars are displayed as in Fig. ^. 
These errors are so small, that their points are almost indistinguishable from the curves. 
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Figure 4: Degeneracy g{E) for the SC 10 x 10 x 10 XY-model obtained from the BHMC- 
fj,G and the BHMC-M methods. Both calculations were performed with AEfix — 2.45. 
The error bars are displayed as in Fig. ^ The curves differ only in a constant value, 
namely ln(5(-0.90)/g(-0.94)). 
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Figure 5: Averaged energy for the SC 10 x 10 x 10 XY-model obtained from the BHMC- 
fiC method (circles), the BHMC-M method (smah squares), and MetropoUs simulations 
(diamonds). Both BHMC calculations were performed with AEfix — 2.45. The values 
obtained from the three procedures are almost indistinguishable. 
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Figure 6: Averaged magnetization for the SC 10 x 10 x 10 XY-model obtained from 
the BHMC-/iC method (circles), the BHMC-M method (small squares) and Metropolis 
simulations (diamonds). Both BHMC calculations were performed with AEfix — 2.45. 
As in Fig. ^, the three curves are almost indistinguishable. 
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and, therefore, it is difficult to obtain estimations for < Ndn{E) > differ- 
ent from zero at tliis energy values. This difficulty increases when AEfi^ 
increases. In the case of the BHMC-M method we have not sampled the 
energy range below E = —0.9 at all. It can be observed that this fact only 
affects the results at low temperatures. For all methods, the errors bars 
become larger around the critical temperature T^, as it would be expected. 

For all these curves, we present error bars of 1.0 standard deviations. For 
eight runs, the confidence level of 99% corresponds to 3.5 standard deviations. 

Since all sets of data obtained from the BHMC-;uC method and the 
Metropolis simulations present error bars of similar size using the same cor- 
relation between samples, we are allowed to compare the speeds of the two 
algorithms. The Metropolis method took 28.5 seconds on a Digital Alpha 
workstation to make one simulation at one temperature point. The BHMC- 
fiC method took 440 seconds to perform one calculation on the whole tem- 
perature range using the same machine. It means that our method calculates 
the whole temperature range 0.7 < T < 4.7 using approximately 15.5 times 
the computer effort required by a Metropolis simulation for one point. In 
addition, the BHMC-M method took 409 seconds per run. 

Figs. ^ and 10 compare the heat capacity obtained from the BHMC-/iC 
and the HMC methods using Metropolis simulations as reference data. Both 
figures display data obtained from the BHMC-/iC and the HMC methods 
using 500 divisions of the energy axis, i.e. AE = 6.0. Fig. ^ presents the 
usual error bars of 1.0 standard deviations and Fig. |10| presents error bars of 
3.5 standard deviations corresponding to a confidence level of 99% for eight 
runs, i. e. the probability of finding the true value outside these error bars 
due only to statistical variations is equal to 1%. It can be observed that all 
methods coincide in the temperature range between T = 2.1 and T = 2.2. 
This is the temperature range of validity of the HMC method according to the 
criteria that the mean energy value calculated at temperature T should be 
less than cr^ away of the mean energy value at Tq, where Tq is the temperature 
used to generate the samples and is the standard deviation on the energy 
of these samples [Q. For temperatures out of this range, the HMC method 
loses its precision but the BHMC-/iC method is still precise on the whole 
temperature range 1.2 < T < 4.7. Since the error bars in the range 2.1 < 
T < 2.2 are similar in size and the correlations between samples are also 
similar, we can compare the speeds of both methods. The HMC method 
took 188 seconds per run on a Digital Alpha workstation. The BHMC-yuC 
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Figure 7: Specific heat for the SC 10 x 10 x 10 XY-model obtained from the BHMC- 
fiC method (circles), the BHMC-M method (smah squares) and Metropohs simulations 
(diamonds). Both BHMC calculations were performed with AEfi^ — 2.45. Within the 
error bars, the curves fit very well onto each other over the temperature range 0.7 < T < 
4.7. The largest differences appear near Tc, as would be expected. This figure gives an 
approximate critical temperature of Tc — 2.16 for the finite size 10 x 10 x 10 system. 




Figure 8: Magnetic susceptibility for the SC 10 x 10 x 10 XY-model obtained from 
the BHMC-/iC method (circles), the BHMC-M method (small squares) and Metropolis 
simulations (diamonds). Both BHMC calculations were performed with AEfix — 2.45. 
As in Fig. ^, the values obtained from the two methods fit very well onto each other, 
showing the largest differences at temperatures near Tc — 2.16. 
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Figure 9: Specific heat for the SC 10 x 10 x 10 XY-model obtained from the BHMC- 
fiC niethod with AEfix — 6.0 (circles) and Metropohs simulations (filled diamonds). The 
curves fit very well onto each other over the temperature range 1.2 <T < 4.7. 
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Figure 10: Specific heat for the SC 10 x 10 x 10 XY-model obtained from the BHMC-/iC 
method with AEfi^ — 6.0 (circles), the HMC method with AE = 6.0 (small squares) and 
Metropolis simulations with 10000 samples (filled diamonds), all methods coincide in the 
temperature range 2.1 < T < 2.2, but the HMC method deviates for temperatures outside 
this range. Instead, the BHMC-/iC method remains precise up to the range 1.2 < T < 4.7 
(see Fig. ^) and using only 2.2 times more computer effort than the HMC method. This 
figure shows error bars of 3.5 standard deviations, i.e. a confidence level of 99% for eight 
runs. 
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method took 411 seconds per run using the same machine. It means that 
our method gives precise results for the temperature range 1.2 < T < 4.7, 
i.e. 35 times larger than the range 2.1 < T < 2.2 of precise results given by 
the HMC method, using only 2.2 times more CPU time. 

7 Conclusions 

We have proposed a way of applying the BHMC method to continuous sys- 
tems. We have found results for the three-dimensional SC XY-model in 
excellent agreement with the HMC method and Metropohs simulations. Our 
method calculates the temperature range 0.7 < T < 4.7 using less computer 
time than the Metropolis simulations for 16 points, and calculates with preci- 
sion the temperature range 1.2 < T < 4.7 using only 2.2 times more computer 
effort than the HMC method for the range 2.1 < T < 2.2. 

The strategy proposed could also be applied to other thermodynamic 
systems with continuous degrees of freedom, for example, the Heisenberg 
model and fluid models. Our strategy can be summarized as follows: 

• Choose a protocol of random movements in the space of states of the 
system, such that for each allowed movement the probability to perform 
it equals the probability to revert it. These movements are only virtual, 
in the sense that they are not executed. They are introduced only to 
calculate Nup, Ndn and, therefore, to estimate the density of states 
g(E). The protocol would change the system from an old io a new state. 
Express the protocol giving the probability density function (p.d.f.) for 
the values of the system variables at the new state. These new values 
are, therefore, random variables. 

• Find the p.d.f. /ab(A£'), i.e. the probability of obtaining an energy 
change between AE and AE + dAE using the chosen protocol. This 
can be done as follows. 

1. Express the energy of the new state Ency,, as a function of the 
new values of the system variables. Enew is, therefore, a random 
variable. 

2. Use the usual rules of finding the p.d.f. for a function of random 
variables to obtain /E„e„(-E), i.e. the p.d.f. of Enew 
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3. Replace E in fE„^^{E) by AE + Eou to obtain, fAE{AE). Here 
Eoid denotes the energy of the old state. 

• Take iV„p {Ndn) for the old state as proportional to f/^E^AEfix) 

ifAE{—^Efix)), and continue with the BHMC method as usual. AEfi^ 
denotes the fixed energy interval of the BHMC method. 

The BHMC method calculates the degeneracy g{E) either by adding 
terms of the form In g{E + AEfi^) — In g{E) or by integrating (3{E), start- 
ing in both cases from a given initial point {Eo,\ng{Eo)). If \ng[Eo) = 
is taken, g{E) can be reinterpreted as the ratio of the probability to obtain 
at random an energy between E and E + dE to the probability to obtain 
an energy between E^ and Eq + dE. This interpretation of the values g{E) 
obtained from the BHMC method is valid both for discrete and continuous 
systems. 

Many different strategies can be used to take the samples in the BHMC 
method. In the present paper we used two of them, namely, a micro-canonical 
sampling process that maintains the system inside a narrow energy window 
(BHMC-/iC) p, and a random walk on the energy axis using Metropolis [0 
steps (BHMC-M) 0. Between these two methods, we prefer the first one, be- 
cause of the following reasons. First, it resembles better the micro-canonical 
character of the averages < Nup{E) > and < Ndn{,E) >. Second, it offers 
a fine control on the correlations between samples and the relaxation time 
at each energy value. Third, two micro-canonical simulations are enough to 
estimate (3{E) at any desired energy and this fact is independent of the sys- 
tem size. Therefore, it seems that the computer time will grow linearly with 
the system size but, of course, the possibility of using data of one simulation 
for two energy points discussed in Sec. || will be lost. This will be studied 
in a future paper. Fourth, many other kinds of micro-canonical simulations 
can be employed to take the samples, different from the accepting-rejecting 
process used in this paper but maintaining also a detailed balance condition 
with equal probabilities for all configurations inside a window. For instance, 
the new configurations can be produced to be always inside the window, and 



this will decrease relaxation and correlation times [|T8], |T9[. Another option 



would be to use Creutz's demons to take the sampleslT 



Other possible sampling strategy is the Entropic Sampling of Ref. P0| , ^ 



To use it in conjunction with the BHMC method is a very interesting idea, 
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because the BHMC method gives direct estimations of lng{E + AEfix) — 
In g{E) and, therefore, no previous estimation of the entropy is needed to 
perform the samphng process. Something similar has been employed in the 
original work of Ref. 0], but estimating these differences only from the last 
sample. 

In conclusion, the strategy proposed here in order to apply the BHMC 
method to systems with continuous degrees of freedom gives excellent results 
for the three dimensional XY-model and seems feasible of being applied to 
other thermodynamic systems. It has to be noted, that the BHMC method 
seems a very precise method to find f3{E) and In g{E) for every thermo- 
dynamic system and it seems also extensible to the calculation of g{E) for 
non-thermodynamic systems. 
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